Language/OS - Multiplatform Resource Library/LANGUAGE OS.iso/a_utils/expanded.lha/expanded/src/MultiMap.C

C/C++ Source or Header | 1992-03-19 | 44KB | 1,638 lines

//
// Linear-Affine-Projective Geometry Package
//
// MultiMap.C
//
// $Header$
//
// William J.R. Longabaugh
// University of Washington
//
// Some routines are courtesy of Steve Mann.
//
// Implementation of the linear-affine-projective geometry
// package described in William J.R. Longabaugh, "An Expanded
// System for Coordinate-Free Geometric Programming", Master's
// thesis, University of Washington, 1992.
//
// Copyright (c) 1992, William J.R. Longabaugh
//   Copying, use and development for non-commercial purposes permitted.
//   All rights for commercial use reserved.
//   This software is unsupported and without warranty; it is
//   provided "as is".
//
// ***********************************************************************

#include "Lap.h"

static int NumSimp(int d, int n);

// *********************************************************************
//
// Returns the value ((d + n) choose n).
// C code provided by Steve Mann.
//

static int NumSimp(int d, int n)
{
  if ( n == -1 ) {
    return 0;
  } else if ( n == 0 ) {
    return 1;
  } else if ( d == 0 ) {
    return 1;
  } else {
    return NumSimp(d - 1, n) + NumSimp(d, n - 1);
  }
}

// *********************************************************************
//
// Multi-indices and tuples of multi-indices are implemented as classes.
//

// Multimaps are currently restricted in degree and dimension.
// These bounds are large enough that they should not create
// difficulties.  The bounds are treated as global constants.

const int DIM_MAX = 10;

class MultiIndex {

  private:

        int size;
        int indices[DIM_MAX];
        Boolean antisym;

  public:

        MultiIndex(void) {};
        MultiIndex(int deg, int dim, Boolean isAnti);
        int SumIndex(void);
        int MItoI(void);
        Boolean NextIndex(void);
    int& Size(void)    {return (size);}
    int& Index(int n)  {return (indices[n]);}
};

// Note the use of SYM_MAX, which is defined in Geom.h:

class MITuple {

  private:

        int size;
        MultiIndex entry[SYM_MAX];
        Boolean antisym;

  public:

    MITuple(MMData& rdat);
    int MITupletoI(void);
    Boolean NextTuple(void);
    int& Size(void)          {return (size);}
    MultiIndex& Entry(int n) {return (entry[n]);}
};

// ***********************************************************************
// ***********************************************************************
//
// MultiIndex Class
//
// ***********************************************************************
// ***********************************************************************
//
// Builds an initial index of degree 'deg' and dimension 'dim'.
// Index is tagged as either symmetric or antisymmetric. Code
// for symmetric case is the C++ version of C code provided
// by Steve Mann.
//

MultiIndex::MultiIndex(int deg, int dim, Boolean isAnti)
{
  antisym = isAnti;
  size = dim - 1;

  if (isAnti) {
    int ones = deg;
    for (int i = size; i >= 0; i--) {
      indices[i] = (ones-- > 0) ? 1 : 0;
    }
  } else {
    indices[size] = deg;
    for (int i = 0; i < size; i++) {
      indices[i] = 0;
    }
  }
}

// *********************************************************************
//
// Returns the sum of all the indices in a multi-index.
// This is a C++ version of C code provided by Steve Mann.
//

int MultiIndex::SumIndex(void)
{
  int sum = 0;

  for (int i = 0; i <= size; i++) {
    sum += indices[i];
  }
  return sum;
}

// *********************************************************************
//
// Determines the associated offset into a linear array for a
// multi-index.  Code for symmetric case is the C++ version of
// C code provided by Steve Mann.
//

int MultiIndex::MItoI(void)
{
  int retval;

  if (antisym) {
    int offset = 0;
    int sum = this->SumIndex();     // degree
    int i = size + 1;               // dimension (remaining slots)
    int readslot = 0;

    while (sum > 0) {
      if (indices[readslot++] == 1) {
        if (i > sum) {
          offset += NumSimp((i - 1 - sum), sum);  // (i - 1 choose sum)
        }
        sum--;
      }
      i--;
    }
    retval = offset;

  } else {
    int sum = this->SumIndex();
    int value = 0;

    for (int i = size; i > 0; i--) {
      sum -= indices[i];
      value += NumSimp(i, sum - 1);
    }
    retval = value;
  }

  return (retval);
}

// *********************************************************************
//
// Transform this multi-index into the next index.  Return TRUE if this
// was done successfully; return FALSE if index is a maximum index.
// Symmetric case code is C++ version of C code provided by Steve Mann.
// No guarantees are made about the state of the index if FALSE is
// returned.
//

Boolean MultiIndex::NextIndex(void)
{
  if (antisym) {
    Boolean can_shift = FALSE;
    int num_shiftable = 0;

//  Start at the most significant slot and cruise down, looking
//  for a 1 that can be shifted up.  Count the number of 1s we
//  pass and change these 1s to zeroes.  When the shift is enabled
//  and a zero is found, do it and then change the appropriate
//  number of zeros back to 1's:

    for (int i = size; i >= 0; i--) {
      if (indices[i] == 1) {
        can_shift = TRUE;
        num_shiftable++;
        indices[i] = 0;
      } else if (can_shift) {
        indices[i] = 1;
        num_shiftable--;
        for (int j = size; num_shiftable > 0; j--) {
          indices[j] = 1;
          num_shiftable--;
        }
        return (TRUE);
      }
    }
    return (FALSE);

  } else {
    int i;

    for (i = size; i >= 0; i--) {
      if (indices[i] != 0) {
        break;
      }
    }
    if (i <= 0) {
      return (FALSE);
    }
    if (i == size) {
      indices[size]--;
      indices[size - 1]++;
    } else {
      indices[size] = indices[i] - 1;
      indices[i] = 0;
      indices[i - 1]++;
    }
    return (TRUE);
  }
}

// ***********************************************************************
// ***********************************************************************
//
// MITuple Class
//
// ***********************************************************************
// *********************************************************************
//
// Creates an initial tuple for given dimension/degree characteristics
// encoded in the MMData object.
//

MITuple::MITuple(MMData& rdat)
{
  antisym = rdat.IsAntisym();
  size = rdat.Numsym();
  for (int i = 0; i < size; i++) {
    entry[i] = MultiIndex(rdat.Deg(i), rdat.Dim(i), antisym);
  }
}

// *********************************************************************
//
// Determine the associated offset into linear storage for this
// multi-index tuple.  The first tuple has index 0.
//

int MITuple::MITupletoI(void)
{
  int value = 0;

// With antisymmetric maps, the implementation is restricted to
// a single symmetry set (i.e. one multi-index):

  if (antisym) {
    value = entry[0].MItoI();
  } else {

//  Start with the rightmost entry in the tuple, and figure out
//  its offset.  Multiply this by the block size.  Do this for
//  each tuple entry, and sum the results.

    int blocksize = 1;
    for (int i = size - 1; i >= 0; i--) {
      value += (entry[i].MItoI() * blocksize);
      blocksize *= NumSimp(entry[i].Size(), entry[i].SumIndex());
    }
  }
  return value;
}

// *********************************************************************
//
//  Given a tuple, make it into the next tuple.  Return TRUE if this
//  was done successfully; return FALSE if tuple is max tuple.
//

Boolean MITuple::NextTuple(void)
{
// Call NextIndex on the rightmost tuple entry.  If we get TRUE,
// we are happy.  If we get FALSE, set the overflowing entry to the
// first index and increment the next entry.

  for (int i = size - 1; i >= 0; i--) {
    if (entry[i].NextIndex()) {
      return TRUE;
    } else {
      int sum = entry[i].SumIndex();
      int misize = entry[i].Size();
      entry[i] = MultiIndex(sum, misize + 1, antisym);
    }
  }
  return FALSE;
}


// ***********************************************************************
// ***********************************************************************
//
// MMData Class
//
// ***********************************************************************
// ***********************************************************************
//
// Need to do memberwise initialization, since Matrix class members need
// to do memberwise initialization:
//

MMData::MMData(MMData& v) : t(v.t)
{
  numsym = v.numsym;
  isAntisym = v.isAntisym;
  for (int k = 0; k < numsym; k++) {
    deg[k] = v.deg[k];
    dim[k] = v.dim[k];
  }
}

// ***********************************************************************
//
// Need to do memberwise assignment, since Matrix class members need
// to do memberwise assignment:
//

MMData& MMData::operator=(MMData& v)
{
  t = v.t;
  numsym = v.numsym;
  isAntisym = v.isAntisym;
  for (int k = 0; k < numsym; k++) {
    deg[k] = v.deg[k];
    dim[k] = v.dim[k];
  }
  return *this;
}

// ***********************************************************************
//
// Output for debugging
//

void MMData::debug_out(ostream &c, int indent)
{
  char *ibloc = new char[indent + 1];
  for (int i = 0; i < indent; i++) {
    *(ibloc + i) = ' ';
  }
  *(ibloc + indent) = '\0';

  c << ibloc << ast;
  c << ibloc << "MMData object\n";
  if (isAntisym) {
    c << ibloc << "Map is antisymmetric\n";
  } else {
    c << ibloc << "Map is partially (or fully) symmetric\n";
  }

  c << ibloc << "Matrix representation of control net:\n";
  t.debug_out(c, indent + ERR_IND);
  c << ibloc << "Number of symmetry sets: " << numsym << "\n";

  c << ibloc << "Degree of each symmetry set:";
  for (i = 0; i < numsym; i++) {
    if (i % 4 == 0) {c << "\n" << ibloc;}
    c << deg[i] << "  ";
  }
  c << "\n";

  c << ibloc << "Dimension of each symmetry set:";
  for (i = 0; i < numsym; i++) {
    if (i % 4 == 0) {c << "\n" << ibloc;}
    c << dim[i] << "  ";
  }
  c << "\n";
  c << ibloc << ast;

  delete ibloc;
  return;
}

// ***********************************************************************
//
// Returns the size of the net needed to store the map:
//

int MMData::NetSize(void)
{
  int total;

// Antisym: (dim choose deg)

  if (isAntisym) {
    total = NumSimp(dim[0] - deg[0], deg[0]);
  } else {

//  For each symmetry set, determine the number of entries.
//  (dim + deg - 1 choose deg)  The total number of slots is
//  the product of these entries:

    total = 1;
    for (int i = 0; i < numsym; i++) {
      total *= NumSimp(dim[i] - 1, deg[i]);
    }
  }
  return total;
}

// ***********************************************************************
//
// Returns the number of the symmetry set corresponding to a particular
// argument nember:
//

int MMData::WhichSet(int argnum)
{
  int inset;
  int symsum = deg[0];
  for (inset = 0; inset < numsym; inset++) {
    if (argnum < symsum) {
      break;
    }
    symsum += deg[inset + 1];
  }
  return (inset);
}

// ***********************************************************************
//
// Updates the MMData to remove the specified symmetry set if it has been
// fully evaluated:
//

void MMData::RemoveSymSet(int inset)
{
  if (deg[inset] == 0) {
    numsym--;
    for (int i = inset; i < numsym; i++) {
      deg[i] = deg[i + 1];
      dim[i] = dim[i + 1];
    }
  }
}

// ***********************************************************************
// ***********************************************************************
//
// MultiMap Class
//
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// Given a multi-index tuple, build an argument of standard basis elements
// for the tuple:
//

void MultiMap::BasisArg(MITuple& tup, Space& s, GeObList* gl)
{
// Step thru each symmetry set.  The value of each field of the
// MultiIndex in that set tells how many times to include each
// standard basis member:

  int slot = 0;
  *gl = GeObList();
  for (int k = 0; k < tup.Size(); k++) {
    Basis curbas = s[slot].StdBasis();
    for (int j = tup.Entry(k).Size(); j >= 0; j--) {
      int n = tup.Entry(k).Index(j);
      for (int m = 0; m < n; m++) {
        *gl = *gl + GeObList(curbas[tup.Entry(k).Size() - j]);
        slot++;
      }
    }
  }
  return;
}

// ***********************************************************************
//
// Get the matrix representation of the nth standard basis object
// (applies to multimaps that are in fact simple maps):
//

Matrix MultiMap::GetStdImage(int n)
{
  static char* sig = "Matrix MultiMap::GetStdImage(int)";

// Make sure the map is in fact a simple map:

  if ((data.Numsym() != 1) || (data.Deg(0) != 1)) {
    errh.ErrorExit(sig, "Multimap not a simple map", *this);
  }

// The matrix of the map encodes the images of the objects in
// the standard basis.  Get the nth image.

  MITuple tup(data);
  for (int k = 0; k < n; k++) {
    tup.NextTuple();
  }

  Matrix retval = data.Net()[tup.MITupletoI()];
  return (retval);
}

// ***********************************************************************
//
// Does a single layer of multimap evaluation for both symmetric and
// antisymmetric cases:
//

void MultiMap::Eval(int argnum, ScalarList& coords, MMData* rdat)
{
  if (rdat->IsAntisym()) {

    // Calculate the sign to apply to this evaluation.  Depends on
    // whether we are evaluating an odd or even argument:

    Scalar oesign =  ((argnum % 2) == 0) ? 1.0 : -1.0;

    // Antisymmetric maps consist of only 1 symmetry set.
    // Generate a first tuple that has the degree of this
    // symmetry set reduced.  Figure out how many new vectors
    // will be generated in this layer of evaluation and build
    // a new matrix.  Evaluate and copy the new matrix back.

    rdat->Deg(0)--;
    Matrix newmat(rdat->NetSize(), rdat->Net().Columns());
    MITuple tup(*rdat);
    Boolean working = TRUE;
    while (working) {
      this->PtEvalAnti(coords, &tup, rdat, &newmat, oesign);
      working = tup.NextTuple();
    }
    rdat->Net() = newmat;

  } else {  // Symmetric maps

    // Figure out which symmetry set the argument belongs to:

    int inset = rdat->WhichSet(argnum);

    // The net will be collapsed within that symmetry set.
    // Generate a first tuple that has the degree of the appropriate
    // symmetry set reduced and evaluate:

    rdat->Deg(inset)--;

    MITuple tup(*rdat);
    Boolean working = TRUE;
    while (working) {
      this->PtEval(coords, &tup, inset, rdat);
      working = tup.NextTuple();
    }

   // Possibly update the symmetry set information:

    rdat->RemoveSymSet(inset);
  }
}

// ***********************************************************************
//
// Evaluates a new net entry for symmetric maps:
//

void MultiMap::PtEval(ScalarList& coords, MITuple* tup,
              int el, MMData* rdat)
{
  Matrix hold = ZeroMatrix(1, rdat->Net().Columns());
  int n = tup->Entry(el).Size();
  int cslot = n;

  for (int i = 0; i <= n; i++) {
    tup->Entry(el).Index(i)++;
    hold += coords[cslot--] * rdat->Net()[tup->MITupletoI()];
    tup->Entry(el).Index(i)--;
  }
  rdat->Net()[tup->MITupletoI()] = hold[0];
}

// ***********************************************************************
//
// Evaluates a new net entry for antisymmetric maps.  Note that results
// are stored in another matrix (antisymmetric evaluation cannot be
// done in place):
//

void MultiMap::PtEvalAnti(ScalarList& coords, MITuple* tup,
                  MMData* rdat, Matrix* newmat, Scalar oesign)
{
  Scalar altsign;
  int onecount = 0;
  Matrix hold = ZeroMatrix(1, rdat->Net().Columns());
  int n = tup->Entry(0).Size();

  for (int i = n; i >= 0; i--) {
    if (tup->Entry(0).Index(i) == 0) {
      tup->Entry(0).Index(i)++;
      altsign = ((onecount % 2) == 0) ? 1.0 : -1.0;
      hold += oesign * altsign *
          coords[n - i] * rdat->Net()[tup->MITupletoI()];
      tup->Entry(0).Index(i)--;
    } else {
      onecount++;
    }
  }
  (*newmat)[tup->MITupletoI()] = hold[0];
}

// ***********************************************************************
//
// Partial evaluation core routine.
//

void MultiMap::PartialEval(GeObList& t, SpaceList* nd, MultiMap* ret)
{
  static char* sig =
      "void MultiMap::PartialEval(GeObList&, SpaceList*, MultiMap*)";

// We have no business trying to evaluate maps that are fully evaluated:

  if (fulleval) {
    errh.ErrorExit(sig, "Map is already fully evaluated", *this, t);
  }

// Transfer the map data into a holding block:

  MMData partial = data;

// Require that the length of the object list match the component count of
// the domain space:

  int numobj = t.Length();
  if (numobj != domain.CPSpaceSize()) {
    errh.ErrorExit(sig,
          "Number of objects does not match cartesian product space size",
          *this, t);
  }

// If at any time we get a tangent space vector in an argument slot
// that expects an affine point, the range of the map, if it is an
// affine space, changes to a tangent space!

  Boolean SwitchRange = FALSE;

// Run through the list of geometric objects, doing one layer of
// evaluation for each one.

  ScalarList coords;
  int evalarg = 0;
  GeOb temp;
  Space tempsp;

  for (int j = 0; j < numobj; j++) {

    // Skip evaluation if the object is NULL.  We need to record the
    // space that was skipped so we know what the new domain space
    // is when we are done:

    if (t[j].Holds() == NULL_GEOB) {
      *nd = *nd + SpaceList(domain[j]);
      evalarg++;
      continue;
    }

    // Cast into appropriate space.  Tangent space vectors can be used instead
    // of affine points:

    GeObType targtype = domain[j].NativeType();
    temp = t[j];
    if (t[j].CanMapTo(targtype)) {
      temp = temp.MapTo(targtype);
      tempsp = temp.SpaceOf();
    } else if ((targtype == AFF_POINT) && temp.CanMapTo(AFF_VECTOR)) {
      temp = temp.MapTo(AFF_VECTOR);
      tempsp = temp.SpaceOf().GetSpace(AFFINE);
      SwitchRange = (rettype == AFF_POINT);
    } else {
      errh.ErrorExit(sig, "Object cannot be mapped to domain space", *this, t);
    }

    if (tempsp != domain[j]) {
      errh.ErrorExit(sig, "Object cannot be mapped to domain space", *this, t);
    }

    // Determine the coordinates of the argument wrt the standard
    // basis of the component space, and pass this off to the single
    // layer evaluation routine:

    coords = (tempsp.StdBasis())(temp);
    this->Eval(evalarg, coords, &partial);
  }

// After this loop, "partial" holds the new data for the multimap,
// but its net matrix is padded with useless data left over
// from the evaluation.  Copy the data to the new map's MMData, only
// keeping the useful stuff.  If the range has switched, the representation
// of the net needs to be truncated:

  Matrix dehoist;
  if (SwitchRange) dehoist = Transpose(range.HoistTangent());
  ret->range = (SwitchRange) ? range.GetSpace(TANGENT) : range;
  ret->data = partial;
  int rowcount = partial.NetSize();
  int colcount = ret->range.MatrixSize();
  ret->data.Net() = Matrix(rowcount, colcount);
  for (int r = 0; r < rowcount; r++) {
    if (SwitchRange) {
      ret->data.Net()[r] = (partial.Net()[r] * dehoist)[0];
    } else {
      ret->data.Net()[r] = partial.Net()[r];
    }
  }

// Figure out odds and ends that are needed to build the new multimap:

  ret->holds = holds;
  ret->fulleval = (evalarg == 0);
  ret->rettype = ret->range.NativeType();
  return;
}

// ***********************************************************************
//
// Dumps the data for a multimap:
//

void MultiMap::data_out(ostream &c, int indent, char* name)
{
  char *ibloc = new char[indent + 1];
  for (int i = 0; i < indent; i++) {
    *(ibloc + i) = ' ';
  }
  *(ibloc + indent) = '\0';

  c << ibloc << ast;
  c << ibloc << name;
  c << ibloc << "Type of multimap: ";
  MultiTypeOut(c, type);
  c << "\n";
  c << ibloc << "Currently holds: ";
  MultiTypeOut(c, holds);
  c << "\n";
  if (holds != NULL_MULTI) {
    c << ibloc << "Domain space:\n";
    domain.debug_out(c, indent + ERR_IND);
    c << ibloc << "Range space:\n";
    range.debug_out(c, indent + ERR_IND);
    c << ibloc << "Domain type: ";
    GeObTypeOut(c, domtype);
    c << "\n";
    c << ibloc << "Return type: ";
    GeObTypeOut(c, rettype);
    c << "\n";
    if (fulleval) {
      c << ibloc << "Map has been fully evaluated\n";
    } else {
      c << ibloc << "Map has not been fully evaluated\n";
    }
    data.debug_out(c, indent + ERR_IND);
  }
  c << ibloc << ast;

  delete ibloc;
  return;
}

// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Need to do memberwise initialization, since class members need
// to do memberwise initialization:
//

MultiMap::MultiMap(MultiMap& v) : range(v.range),
                  domain(v.domain), data(v.data)
{
  holds = v.holds;
  type = ANY_MULTI;
  fulleval = v.fulleval;
  domtype = v.domtype;
  rettype = v.rettype;
}

// ***********************************************************************
//
// Need to do memberwise assignment, since class members need
// to do memberwise assignment:
//

MultiMap& MultiMap::operator=(MultiMap& v)
{
//
// This weird checking is required to bypass the apparent inheritance of
// assignment under g++ 1.37:
//
  if ((type != ANY_MULTI) && (type != v.holds) && (v.holds != NULL_MULTI)) {
    errh.ErrorExit("MultiMap& MultiMap::operator=(MultiMap&)",
                   "Illegal assignment attempted", *this, v);
  }
  range = v.range;
  domain = v.domain;
  data = v.data;
  holds = v.holds;
  fulleval = v.fulleval;
  domtype = v.domtype;
  rettype = v.rettype;
  return *this;
}

// ***********************************************************************
//
// Output for debugging:
//

void MultiMap::debug_out(ostream &c, int indent)
{
  static char* name = "MultiMap Object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
// Converts a linear or affine map into a multimap:
//

MultiMap::MultiMap(Map &s)
{
  static char* sig = "MultiMap::MultiMap(Map&)";

  type = ANY_MULTI;

// The map must be linear or an affine map from an entire affine space:

  if ((s.Holds() != LIN_MAP) && (s.Holds() != AFF_MAP)) {
    errh.ErrorExit(sig, "Map must be linear or affine", s);
  }

  if ((s.Holds() == AFF_MAP) &&
      (s.Domain().EmbeddingSpace().Holds() != AFF_SPACE)) {
    errh.ErrorExit(sig, "Domain of affine map must be an affine space", s);
  }

// The map must be from a full space:

  if (!s.Domain().IsFullSpace()) {
    errh.ErrorExit(sig, "Domain must be an entire space", s);
  }

  range = s.Range().EmbeddingSpace();
  domain = s.Domain().EmbeddingSpace();

  if ((range.Holds() != VEC_SPACE) && (range.Holds() != AFF_SPACE)) {
    errh.ErrorExit(sig, "Range must be an affine or vector space", range);
  }

// The domain must be an atomic space:

  if (domain.CPSpaceSize() != 1) {
    errh.ErrorExit(sig, "Domain must be an atomic space", s);
  }

  holds = (s.Holds() == AFF_MAP) ? MULTI_AFFINE : MULTI_LINEAR;
  fulleval = (domain.Holds() == NULL_SPACE);
  domtype = s.DomainType();
  rettype = s.ReturnType();
  data.IsAntisym() = FALSE;

// These maps have only one symmetry set:

  data.Numsym() = 1;
  data.Deg(0) = 1;
  data.Dim(0) = domain.MatrixSize();

  if (data.Dim(0) > DIM_MAX) {
    errh.ErrorExit(sig, "Multimap dimension restriction exceeded", s);
  }

// The matrix of the map encodes the images of the objects in
// the standard basis.

  Matrix temp = s.MatrixRep();
  MITuple tup(data);
  data.Net() = Matrix(data.Dim(0), temp.Columns());

  for (int k = 0; k < data.Dim(0); k++) {
    data.Net()[tup.MITupletoI()] = temp[k];
    tup.NextTuple();
  }
}

// ***********************************************************************
//
// Full evaluation.
//

GeOb MultiMap::operator()(GeOb& v)
{
  static char* sig = "GeOb MultiMap::operator()(GeOb&)";

// We have no business trying to evaluate maps that are fully evaluated:

  if (fulleval) {
    errh.ErrorExit(sig, "Map is already fully evaluated", *this, v);
  }

// Cast the argument into the domain space.  If the argument is a
// vector tuple and we expect a point tuple, can still go:

  GeOb temp;
  Space tempsp;
  GeObType ret;
  Space rng;
  Boolean dehoist = FALSE;

  if (v.CanMapTo(domtype)) {
    temp = v.MapTo(domtype);
    tempsp = temp.SpaceOf();
    ret = rettype;
    rng = range;
  } else if ((domtype == AFF_POINT) && v.CanMapTo(AFF_VECTOR)) {
    temp = v.MapTo(AFF_VECTOR);
    tempsp = temp.SpaceOf().GetSpace(AFFINE);
    if (rettype == AFF_POINT) {
      ret = AFF_VECTOR;
      rng = range.GetSpace(TANGENT);
      dehoist = TRUE;
    } else {
      ret = rettype;
      rng = range;
    }
  } else {
    errh.ErrorExit(sig, "Object cannot be mapped to domain space", *this, v);
  }

  if (tempsp != domain) {
    errh.ErrorExit(sig, "Object cannot be mapped to domain space", *this, v);
  }

// Step through each component of the argument.  For each component,
// determine the coordinates the the argument wrt the standard
// basis of the component space, and pass this off to the single
// layer evaluation routine:

  int n = tempsp.CPSpaceSize();
  MMData partial = data;

  for (int j = 0; j < n; j++) {
    ScalarList coords = (tempsp[j].StdBasis())(temp[j]);
    this->Eval(0, coords, &partial);
  }

// After this loop, "partial" holds the std. coords of the object
// in the range space.  Build the GeOb from this information.
// We may need to drop the representation down to the tangent space:

  Matrix retmat(partial.Net()[0]);
  if (dehoist) {
    retmat *= Transpose(range.HoistTangent());
  }
  GeOb retval(ret, rng, retmat);
  return retval;
}

// ***********************************************************************
//
// Partial evaluation, where we need to build a new CP Space for the
// new domain:
//

MultiMap MultiMap::operator()(GeObList &t)
{
  static char* sig = "MultiMap MultiMap::operator()(GeObList&)";
  MultiMap retval;
  SpaceList ndlist;

// Call the core routine to do the grunt work:

  this->PartialEval(t, &ndlist, &retval);

// Build the required domain space:

  if (retval.fulleval) {
    retval.domain = Space();
  } else if (ndlist.Length() == 1) {
    retval.domain = ndlist[0];
  } else if (holds == MULTI_LINEAR) {
    retval.domain = VSpace("Domain for partially evaluated MLM",
               ndlist, FALSE);
  } else {
    retval.domain = ASpace("Domain for partially evaluated MAM",
               ndlist, FALSE);
  }
  retval.domtype = retval.domain.NativeType();

  return retval;
}

// ***********************************************************************
//
// Partial evaluation with specified target domain space:
//

MultiMap MultiMap::operator()(Space& newdom, GeObList &t)
{
  static char* sig = "MultiMap MultiMap::operator()(Space&, GeObList&)";
  MultiMap retval;
  SpaceList ndlist;

// Call the core routine to do the grunt work:

  this->PartialEval(t, &ndlist, &retval);

// Make sure the new domain space provided matches the one needed:

  int num = ndlist.Length();
  if (num != newdom.CPSpaceSize()) {
      errh.ErrorExit(sig,
             "Target domain space size mismatch", *this, newdom, t);
    }
  for (int j = 0; j < num; j++) {
    if (ndlist[j] != newdom[j]) {
      errh.ErrorExit(sig, "Target domain space mismatch", *this, newdom, t);
    }
  }

// Finish building the new map:

  retval.domain = newdom;
  retval.domtype = newdom.NativeType();

  return retval;
}

// ***********************************************************************
// ***********************************************************************
//
// MLM Class
//
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// Creates the standard multi-linear maps for Euclidean spaces.
//

MLM::MLM(StdMaps mpt, VSpace& sp)
{
  type = MULTI_LINEAR;
  switch (mpt) {
    case INNER_PROD:
      this->DotProduct(sp);
      break;
    case CROSS_PROD:
      this->CrossProduct(sp);
      break;
    default:
      errh.ErrorExit("MLM::MLM(StdMaps, VSpace&)", "Unexpected error");
      break;
  }
}

// ***********************************************************************

void MLM::DotProduct(VSpace& sp)
{
  static char* sig = "void MLM::DotProduct(VSpace&)";

// Build data info.  There is currently a size restriction on MultiMaps:

  if (sp.Dim() > DIM_MAX) {
    errh.ErrorExit(sig, "Multimap dimension restriction exceeded", sp);
  }
  data.Numsym() = 1;
  data.IsAntisym() = FALSE;
  data.Deg(0) = 2;
  data.Dim(0) = sp.Dim();

// Fill in the data matrix.  Entries are either 0 or 1 (vectors in the
// space Reals).  The slot in the data matrix depends on the multi-index
// ordering:

  MITuple tup(data);
  int count = data.NetSize();
  data.Net() = Matrix(count, 1);

  for (int k = 0; k < count; k++) {
    data.Net()[tup.MITupletoI()] = this->DotResult(tup);
    tup.NextTuple();
  }

  holds = MULTI_LINEAR;
  fulleval = FALSE;
  domain = VSpace("Argument space for dot product", SpaceList(sp, sp), FALSE);
  range = Reals;
  domtype = domain.NativeType();
  rettype = Reals.NativeType();
}

// ***********************************************************************
//
// Returns the image corresponding to the given tuple for the standard
// Euclidean dot product.
//

Scalar MLM::DotResult(MITuple& tup)
{
// The tuple is only going to have one symmetry group.  The maximum
// degree is 2.  Scan the MultiIndex.  If we hit a 2, the arguments
// match and we need a 1.0.  Otherwise, the result is 0.0:

  int i;
  int top = tup.Entry(0).Size();

  for (i = 0; i <= top; i++ ) {
    int t = tup.Entry(0).Index(i);
    if (t == 1) {
      return (0.0);
    } else if (t == 2) {
      return (1.0);
    }
  }
}

// ***********************************************************************
//
// For n-dimensional space V, cross product is (V x V x ... x V) -> V,
// where the domain space repeats V (n - 1) times.
//

void MLM::CrossProduct(VSpace& sp)
{
  static char* sig = "void MLM::CrossProduct(VSpace&)";

// Build data info.  There is currently a size restriction on MultiMaps:

  if (sp.Dim() > DIM_MAX) {
    errh.ErrorExit(sig, "Multimap dimension restriction exceeded", sp);
  }
  data.Numsym() = 1;
  data.IsAntisym() = TRUE;
  data.Deg(0) = sp.Dim() - 1;
  data.Dim(0) = sp.Dim();

// Fill in the data matrix.  Entries are either plus or minus the
// basis vectors in V.  The slot in the data matrix depends on the multi-index
// ordering:

  MITuple tup(data);
  int count = data.NetSize();
  data.Net() = Matrix(count, sp.MatrixSize());

  for (int k = 0; k < count; k++) {
    // The following line produces a g++ error, so store result in a temp:
    // data.Net()[tup.MITupletoI()]=(this->CrossResult(tup, sp.StdBasis()))[0];
    Matrix holdmat = this->CrossResult(tup, sp.StdBasis());
    data.Net()[tup.MITupletoI()] = holdmat[0];
    tup.NextTuple();
  }

  holds = MULTI_LINEAR;
  fulleval = (sp.Dim() == 1);
  if (fulleval) {
    domain = Space();
  } else if (sp.Dim() == 2) {
    domain = sp;
  } else {
  domain = VSpace("Argument space for cross product",
          SpaceList(sp, sp.Dim() - 1), FALSE);
  }
  range = sp;
  domtype = domain.NativeType();
  rettype = range.NativeType();
}

// ***********************************************************************
//
// Returns the image corresponding to the given tuple for the standard
// right-handed Euclidean cross product.
//

Matrix MLM::CrossResult(MITuple& tup, Basis& stdbas)
{
// The tuple is only going to have one symmetry group, with one
// zero.  Find the slot with the zero, get the corresponding basis
// vector.  The sign is determined by the number of swaps needed to
// get the monotonic ordering, e.g. 12345:

  int top = tup.Entry(0).Size();
  Scalar swapsign;

  for (int i = top; i >= 0; i--) {
    int t = tup.Entry(0).Index(i);
    if (t == 0) {
      swapsign = (i % 2 == 0) ? 1.0 : -1.0;
      // G++ 1.37.1 manages to botch this next statement up by apparently
      // destructing the temporary GeOb before copying the matrix:
      // Matrix retmat = stdbas[top - i].MatrixRep() * swapsign;
      // Fix it by assigning the GeOb to a temporary:
      GeOb holdg = stdbas[top - i];
      Matrix retmat = holdg.MatrixRep() * swapsign;
      return (retmat);
    }
  }
}

// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Used to cast GeObs to MLMs.  Too large to inline under g++ 1.37.1:
//

MLM::MLM(GeOb& s)
{
  type = MULTI_LINEAR;
  *this = MLM(MultiMap(s));
}

// ***********************************************************************
//
// Used to cast a general multimap down to a multi-linear map:
//

MLM::MLM(MultiMap& m): (m)
{
  type = MULTI_LINEAR;
  if ((holds != MULTI_LINEAR) && (holds != NULL_MULTI)) {
      errh.ErrorExit("MLM::MLM(MultiMap&)",
          "Attempted to cast a non-multi-linear map to a multi-linear map", m);
  }
}

// ***********************************************************************
//
// Need to do memberwise assignment, since class members need
// to do memberwise assignment:
//

MLM& MLM::operator=(MLM& v)
{
  range = v.range;
  domain = v.domain;
  data = v.data;
  holds = v.holds;
  fulleval = v.fulleval;
  domtype = v.domtype;
  rettype = v.rettype;
  return *this;
}

// ***********************************************************************
//
// Output for debugging:
//

void MLM::debug_out(ostream &c, int indent)
{
  static char* name = "MLM Object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
// Build a multi-linear map:
//

MLM::MLM(VSpace& s1, IntList& symmetry, BasisList& bases,
     VSpace& s2, GeObList& vectors)
{
  static char* sig =
     "MLM::MLM(VSpace&, IntList&, BasisList&, VSpace&, GeObList&)";
  int m;

// Odds & ends:

  type = MULTI_LINEAR;
  holds = MULTI_LINEAR;
  domain = s1;
  range = s2;
  domtype = domain.NativeType();
  rettype = range.NativeType();
  fulleval = (domain.Holds() == NULL_SPACE);

// Check validity of symmetry specification.  There is a symmetry
// set size restriction on multimaps:

  int num = symmetry.Length();
  if (num > SYM_MAX) {
    errh.ErrorExit(sig, "Too many symmetry sets specified", symmetry);
  } else {
    data.Numsym() = num;
  }
  if (num != bases.Length()) {
    errh.ErrorExit(sig, "Incorrect number of bases specified",
           symmetry, bases);
  }

// Go through the symmetry specification entry by entry:

  data.IsAntisym() = FALSE;
  int cpsize = 0;
  for (int i = 0; i < num; i++) {
    int val = symmetry[i];

    // Handle the specification of antisymmetric maps:

    if (val < 0) {
      if (num > 1) {
        errh.ErrorExit(sig, "Antisymmetric map can only have 1 symmetry set",
               symmetry);
      } else {
    data.IsAntisym() = TRUE;
    val = -val;
      }
    } else if (val == 0) {
      errh.ErrorExit(sig, "Symmetry list entries cannot equal 0", symmetry);
    }

    // Check validity within current symmetry set:

    if ((cpsize + val) > s1.CPSpaceSize()) {
      errh.ErrorExit(sig,
    "Sum of symmetry list is more than number of domain space components",
    s1, symmetry);
      }
    Space curspace = s1[cpsize];
    for (int j = 1; j < val; j++) {
      if (s1[cpsize + j] != curspace) {
        errh.ErrorExit(sig,
        "Domain component spaces not identical within a symmetry set",
         s1, symmetry);
      }
    }
    if (bases[i].SpaceOf() != curspace) {
      errh.ErrorExit(sig,
        "Specified basis not in corresponding domain component space",
         s1, curspace, bases, symmetry);
    }
    cpsize += val;
    if (curspace.Dim() > DIM_MAX) {
      errh.ErrorExit(sig, "Multimap dimension restriction exceeded", curspace);
    }
    data.Deg(i) = val;
    data.Dim(i) = curspace.Dim();
  }

// Sum of symmetry spec. list must equal the number of components
// in domain space:

  if (cpsize != s1.CPSpaceSize()) {
    errh.ErrorExit(sig,
     "Sum of symmetry list must equal number of domain space components",
     domain, symmetry);
  }

// Fill in a temporary data set that will hold the representations of
// the image objects with respect to the user-specified bases.
// Need to coerce the image elements into the specified range space.
// The slot in the data matrix depends on the multi-index ordering:

  int vecnum = vectors.Length();
  data.Net() = Matrix(vecnum, range.MatrixSize());
  GeOb temp;
  MMData userspec = data;
  MITuple tup(data);

  if (vecnum != data.NetSize()) {
    errh.ErrorExit(sig, "Incorrect number of image vectors specified",
                s1, symmetry, bases, s2, vectors);
    }
  for (m = 0; m < vecnum; m++) {
    if (vectors[m].CanMapTo(rettype)) {
      temp = vectors[m].MapTo(rettype);
      if (temp.SpaceOf() != range) {
        errh.ErrorExit(sig, "Object cannot be mapped to range space",
               range, vectors[m]);
      }
    } else {
      errh.ErrorExit(sig, "Object cannot be mapped to range space",
             range, vectors[m]);
    }
    userspec.Net()[tup.MITupletoI()] = temp.MatrixRep()[0];
    tup.NextTuple();
  }

// Now build a final data set row by row by calculating the images
// correponding to the standard basis elements:

  GeObList basisarg;
  tup = MITuple(data);
  for (m = 0; m < vecnum; m++) {
    MMData partial = userspec;
    this->BasisArg(tup, domain, &basisarg);
    for (int j = 0; j < basisarg.Length(); j++) {
      ScalarList coords = bases[userspec.WhichSet(j)](basisarg[j]);
      this->Eval(0, coords, &partial);
    }
    data.Net()[tup.MITupletoI()] = partial.Net()[0];
    tup.NextTuple();
  }
}

// ***********************************************************************
// ***********************************************************************
//
// MAM Class
//
// ***********************************************************************
// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Used to cast a general multimap down to an MAM:
//

MAM::MAM(MultiMap& m): (m)
{
  type = MULTI_AFFINE;
  if ((holds != MULTI_AFFINE) && (holds != NULL_MULTI)) {
    errh.ErrorExit("MAM::MAM(MultiMap&)",
          "Attempted to cast a non-multi-affine map to a multi-affine map", m);
  }
}

// ***********************************************************************
//
// Need to do memberwise assignment, since class members need
// to do memberwise assignment:
//

MAM& MAM::operator=(MAM &v)
{
  range = v.range;
  domain = v.domain;
  data = v.data;
  holds = v.holds;
  fulleval = v.fulleval;
  domtype = v.domtype;
  rettype = v.rettype;
  return *this;
}

// ***********************************************************************
//
// Output for debugging:
//

void MAM::debug_out(ostream &c, int indent)
{
  static char* name = "MAM object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
// Build a multi-affine map:
//

MAM::MAM(ASpace& s1, IntList& symmetry, BasisList& simplices,
     Space& s2,  GeObList& images)
{
  static char* sig =
     "MAM::MAM(ASpace&, IntList&, BasisList&, Space&, GeObList&)";
  int m;

// Odds & ends:

  type = MULTI_AFFINE;
  holds = MULTI_AFFINE;
  domain = s1;
  range = s2;
  if ((range.Holds() != VEC_SPACE) && (range.Holds() != AFF_SPACE)) {
    errh.ErrorExit(sig, "Range must be an affine or vector space", range);
  }
  domtype = domain.NativeType();
  rettype = range.NativeType();
  fulleval = (domain.Holds() == NULL_SPACE);
  data.IsAntisym() = FALSE;

// Check validity of symmetry specification.  There is a symmetry
// set size restriction on multimaps:

  int num = symmetry.Length();
  if (num > SYM_MAX) {
    errh.ErrorExit(sig, "Too many symmetry sets specified", symmetry);
  } else {
    data.Numsym() = num;
  }
  if (num != simplices.Length()) {
    errh.ErrorExit(sig, "Incorrect number of bases specified",
           symmetry, simplices);
  }

// Go through the symmetry specification entry by entry:

  int cpsize = 0;
  for (int i = 0; i < num; i++) {
    int val = symmetry[i];
    if (val <= 0) {
      errh.ErrorExit(sig, "Symmetry list entries must be >= 1", symmetry);
    }

    // Check validity within current symmetry set:

    if ((cpsize + val) > s1.CPSpaceSize()) {
      errh.ErrorExit(sig,
    "Sum of symmetry list is more than number of domain space components",
    s1, symmetry);
      }
    Space curspace = s1[cpsize];
    for (int j = 1; j < val; j++) {
      if (s1[cpsize + j] != curspace) {
        errh.ErrorExit(sig,
        "Domain component spaces not identical within a symmetry set",
         s1, symmetry);
      }
    }
    if (simplices[i].SpaceOf() != curspace) {
      errh.ErrorExit(sig,
        "Specified basis not in corresponding domain component space",
         s1, curspace, simplices, symmetry);
    }
    if (simplices[i].Holds() != SIMPLEX) {
      errh.ErrorExit(sig, "Specified basis not a simplex",
             simplices, simplices[i]);
    }
    cpsize += val;
    if (curspace.Dim() >= DIM_MAX) {
      errh.ErrorExit(sig, "Multimap dimension restriction exceeded", curspace);
    }
    data.Deg(i) = val;
    data.Dim(i) = curspace.Dim() + 1;
  }

// Sum of symmetry spec. list must equal the number of components
// in domain space:

  if (cpsize != s1.CPSpaceSize()) {
    errh.ErrorExit(sig,
     "Sum of symmetry list must equal number of domain space components",
     domain, symmetry);
  }

// Fill in a temporary data set that will hold the representations of
// the image objects with respect to the user-specified bases.
// Need to coerce the image elements into the specified range space.
// The slot in the data matrix depends on the multi-index ordering:

  int numimg = images.Length();
  data.Net() = Matrix(numimg, range.MatrixSize());
  GeOb temp;
  MMData userspec = data;
  MITuple tup(data);

  if (numimg != data.NetSize()) {
    errh.ErrorExit(sig, "Incorrect number of image objects specified",
                s1, symmetry, simplices, s2, images);
    }
  for (m = 0; m < numimg; m++) {
    if (images[m].CanMapTo(rettype)) {
      temp = images[m].MapTo(rettype);
      if (temp.SpaceOf() != range) {
        errh.ErrorExit(sig, "Object cannot be mapped to range space",
               range, images[m]);
      }
    } else {
      errh.ErrorExit(sig, "Object cannot be mapped to range space",
             range, images[m]);
    }
    userspec.Net()[tup.MITupletoI()] = temp.MatrixRep()[0];
    tup.NextTuple();
  }

// Now build a final data set row by row by calculating the images
// correponding to the standard basis elements:

  GeObList basisarg;
  tup = MITuple(data);
  for (m = 0; m < numimg; m++) {
    MMData partial = userspec;
    this->BasisArg(tup, domain, &basisarg);
    for (int j = 0; j < basisarg.Length(); j++) {
      ScalarList coords = simplices[userspec.WhichSet(j)](basisarg[j]);
      this->Eval(0, coords, &partial);
    }
    data.Net()[tup.MITupletoI()] = partial.Net()[0];
    tup.NextTuple();
  }
}

// ***********************************************************************